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ABSTRACT 

When thin accretion disks around black holes are perturbed, the main restor- 
ing force is gravity. If gas pressure, magnetic stresses, and radiation pressure 
are neglected, the disk remains thin as long as orbits do not intersect. Inter- 
sections would result in pressure forces which limit the growth of perturbations. 
We find that a discrete set of perturbations is possible for which orbits remain 
non-intersecting for arbitrarily long times. These modes define a discrete set of 
frequencies. We classify all long-lived perturbations for arbitrary potentials and 
show how their mode frequencies are related to pattern speeds computed from 
the azimuthal and epicyclic frequencies. We show that modes are concentrated 
near radii where the pattern speed has vanishing radial derivative. 

We explore these modes around Kerr black holes as a possible explanation 
for the high-frequency quasi-periodic oscillations of black hole binaries such as 
GRO J1655-40. The long-lived modes are shown to coincide with diskoseismic 
waves in the limit of small sound speed. While the waves have long lifetime, they 
have the wrong frequencies to explain the pairs of high-frequency quasi-periodic 
oscillations observed in black hole binaries. 

Subject headings: accretion disks - black hole physics - relativity - X-rays:binaries 
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Introduction 



X-ray emission from accretion disks around compact objects often shows high-frequency 
flickering in the form of quasi-periodic oscillations (QPO). Several galactic binaries containing 
black holes show pairs of QPO of fixe d, relatively high frequencies having ratios very close 
to 3:2 jRemillard fcMcClintocl] bood ) . These include GRO J1655 -40 with 300 and 450 Hz 
QPO Jstrohmaveil boOlh: XTE J 1550-564 w ith 184 and 276 Hz jMiller et aliboOlh : GRS 
1915+105 with 113 and 168 Hz jRemillardI [200I : and H1743-322 with 165 and 240 Hz 
faoman et al1l2005h . 



The origin of these quasi-periodic oscillations is still unclear. One possibility is that 
orbiting blobs or hot spots follow non-circular orbits whose epicyclic frequencies are related 
to the obse rved frequencies (Stella, Vietr i fc M orsin k 1999). Thi s mod el has been investigated 
in detail by lSchnittman fc Bertschingerl (120041 ) and ISchnittmaru (120051 ). who use a ray-tracing 
code to model the power spectrum of light curve variations. The model has two serious 
shortcomings. First, special orbits must be chosen by hand (e.g., the radii of the nearly 
circular orbits) to produce the correct frequencies, so this model lacks predictive power. 
Moreover, the blobs can persist only a few rotation periods before being sheared apart by 
differential rotation. 

Another mo d el tha t can reproduce a 3:2 frequency ratio is the relativistic torus model 
of iRezzoUa et al.l (120031 ). To produce the 3:2 ratio, the model requires a black hole spin close 
to maximal. The actual value of frequency depends on the radial extent of the torus. This 
is a free parameter of the model, which therefore also lacks predictive power. 

The fixed frequencies and 3:2 ratio of these QPO pairs suggests that they may repre- 
sent nonlinearly driven resonant frequencies in the curved spacetime around the black hole 
(Abramowicz & Kluzniak 2001; Rebusco 2004; Lee, Abramowicz Sz Kluzniak 2004). This 
idea has been studied using a mathematical model for parametric resonance but has not yet 
been shown to arise automatically from nonlinear dynamics in general relativity. 

Another possible explanation for the fixed frequencies of high-fr equency QPO is disko- 
seismic waves similar to helioseismic (7- modes ( iWagoner et al.ll200ll ). This model has the 
advantage of predicting a discrete set of frequencies. However, a frequency ratio of 3:2 only 
obtains for a few particular values of the black hole spin that are unlikely to obtain for 
all observed systems. Thus, similarly to the hot spot model, this model requires tuning of 
parameters. 

The models can be classified as to whether the oscillating structures are particle-like 
(blobs, with or without resonant interactions) or wave-like (oscillating tori, diskoseismic 
waves). Blob models are in many respects the simplest, as the only important physics is 
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gravity. It is natural to ask whether a similarly simple description is possible for waves. 

In this paper we show that waves in an accretion disk can also be described using simple 
kinematics in curved spacetime. We investigate nearly kinematic density waves, which are 
waves of compression and rarefaction of particles freely orbiting in a disk. These density 
waves are very similar to the density waves which underlie the Lin-Shu hypothesis of galaxy 
spiral structure, apart from several important differences. Our model neglects self-gravity, 
pressure, magnetic stresses, and all forces except the gravity of the central compact object. 
A kinematic model of perturbation s in thin accretion disks around Kerr black holes was 



considered in lAmin fc Frolovl (120061 ) . However, they allow for fluid streamlines to intersect 



forming shocks, as we will describe below. 

Early models for galactic spiral structure (Lindblad 1957, Kalnajs 1973) failed to pre- 
dict sufficiently long-lived patterns because differential rotation relatively quickly winds up 
material arms; moreover, these models neglected self-gravity which is invalid for galactic 
disks (Toomre 1977; Binney & Tremaine 1987). However, accretion disks around black holes 
contain relatively much less mass than their central "bulge" , so self-gravity is much less im- 
portant. If other forces are unimportant — as might be expected for a thin disk — then the 
remaining challenge is to prevent the waves from becoming tig htly wound in a short time. 



This is exactly the problem considered in the galactic context by lKalnajd (Il973l ). Our model 
is a generalization of his kinematic density waves to arbitrary potentials and to the strongly 
curved spacetime around a Kerr black hole. 

We choose to investigate a thin disk, since it is analytically tractable. In this case 
pressure forces are not important, otherwise the disk would not remain very thin. For 
consistency, shocks must not form in the gas. Thus the perturbations must evolve so that 
fluid streamlines avoid intersections. The boundary between regions where streamlines cross 
and where they do not is a caustic surface. Our model attempts to flnd long-lived kinematic 
density waves free from caustics. 



2. Disk Perturbations in Generic Potentials 

In order to study the dynamics of thin accretion disks, we will perturb an initially sta- 
tionary and axisymmetric thin disk and then impose the condition that particle trajectories 
do not intersect. The surface number density is no(f) where f is the radius of an unper- 
turbed orbit. This assumes that the effective potential is axisymmetric and each orbit has a 
conserved angular momentum as well as energy. 

Perturbations are added using the Lagrangian framework based on individual particle 
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trajectories. Thus, at each instant of time, there exists a nonhnear mapping that takes the 
coordinates of a particle in the unperturbed disk to the coordinates of the same particle in 
a perturbed disk. We choose to map the Lagrangian coordinates (which we choose to call f 
and (p) of the particles in the unperturbed disk to the Eulerian spherical coordinates (r, (p, 6) 
of the same particles in the perturbed disk; 9 = 7t/2 for the unperturbed disk. Thus, the 
map is given by 



where R, P, and T are functions we will parameterize below. 

We require the mapping to be smooth and single- valued so that there are no intersections 
of particle orbits, which would lead to shock waves and invalidate our assumption that the 
disk remains cold and thin. In other words, the perturbed disk is given by a diffeomorphism 

of the unperturbed disk. The number density n(r, 0, 6; t) of the perturbed disk follows from 
the transformation of the disk area element rfS. In the absence of intersecting orbits, the 
number of particles dN contained in each surface clement is constant: dN = no{f)df^ ~ 
n{r, 4>, 9; t)dT^. The surface number density of the perturbed disk is then given by 



Here rfS is an area element in the unperturbed disk. The requirement that no intersections 
occur translates as a condition on liS: ^ 0. 

If the disk perturbations were coplanar displacements in fiat space, the ratio dH/dTj 
would be a simple Jacobian determinant. However, the displacements are (in general) three- 
dimensional and we allow for space curvature. Thus a more careful calculation of the area 
element is needed. We do this using an orthonormal basis constructed on shces of constant 
coordinate time {dt — 0). The line element on each time slice is dl'^ — dx^dx^gij — dx^dx^Sij, 
where the hatted indices give the components of a vector in the orthonormal basis and Latin 
indices are summed over the spatial components (with positive signature) . Here, gij are the 
spatial components of the full spacetime metric. 

Orthonormal and coordinate components of three- vectors are related by a triad e*^-: 

dx^ — e^jdx^. Prom the invariance of dP, we find that e\e-'i6ij = gki- For a diagonal spatial 
metric gij (i.e., an orthogonal coordinate basis) we can orient the triad so that its only non- 
zero components are e\ = (no summation). Throughout the rest of the paper, we will 
assume that the spatial metric is diagonal (as it is for the Kerr metric in Boyer-Lindquist 
coordinates) . 

Let us call the three spatial coordinates (r, 0, 9) where 9 — n/ 2 is the usual equatorial 
surface. The components of the infinitesimal displacement vector v— dx, in the right-handed 



r{t) = R{f, 0; t), <P{t) = P(f, 0; t), 9{t) = r(f, 0; t) , 





(2) 
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coordinate basis {dr, d^, -de) are = (dr, d(j), —dO). Using the triad, the vector components 
in the orthonormal basi^ (f, 0, 6) are therefore given by 

= {Adr, Bd(f), -Cdd) , A = ^Jg^r , B = ^/g^ , C = y/g^ . (3) 

The factors {A, B, C) in general can be functions of (r, (p, 6). Every vector is associated with 
a point (r, 0, 6) in the perturbed disk and (f, 0) in the unperturbed disk. 

Let us construct two displacement vectors X and Y whose origin is the point (f , 0, 6) 
in the plane of the unperturbed disk as follows, 

X = Adfr, Y = Bdfct) , (4) 

where A = A{f,(f),6), and analogously for B. Under the diffeomorphism ([ID, we find that 
the vectors are mapped to 

X = a(^] dfr + B(^) df^-c(^) dfO 

Y = a( ^\ d^r + B( ^] df4>-c( ^\ dfe . (5) 

Note that the origin of the vectors is changed under the map, and therefore A,B,C are no 
longer barred. 

The area of the parallelogram in the unperturbed disk formed by X and Y is given by 
dTj = |X X 1^1 = ABdfdcf). This parallelogram is mapped under ([1]) to a parallelogram with 
area dT, = \X ^Y\. Note the vector X x is normal to the perturbed disk. Thus, we have 



X X r 1 



|X X 1^1 ABdfdcj) 



r 4) e 

A{dR/df)df B{dP/df)df -C{dT/df)df 
A{dR/d^)d^ B{dP/d^)d^ -C{dT/d^)d^ 



BC ^ ^ AC ^ AB ^ ^ 

JpTf + -J^JRTfP + -J^JrpO , (6) 



AB AB AB 

and therefore, from ([2]), 

where Jmn is the Jacobian of any functions M(f, 0) and A^(f, 0). 



^Note that we define to point in the —dg direction 



- 6 - 



In order to avoid caustics, we must impose that {riQ/nY 7^ in the disk at each moment 
in time. For the spacetimes of interest, the triad factors (A, 5, C) are nonzero in the disk, 
implying that at least one of the three differential inequalities Jmn 7^ must always be 
satisfied at each point of the disk. In order to systematically solve these inequalities, we will 
work in the epicyclic approximation. Throughout the rest of the paper we will assume that 
A, B and C do not vanish in the plane of the disk. 



r{t) 


= R{r,^ 


b;t) 


= r — e(r, $) cos[uJr\ 


r)t + xij, 




m 


= P{r,^ 


M) 


= + 0o(^, 0) 


;f)t + /3(f 


)e{f, $) sm[uJr{f)t + xir, $)] 


e{t) 


= T{f,(, 


M) 


= 9 — 7(f, $) cos[ujg 


(f)t + C(f, 


$)] . 



3. Epicyclic Approximation 

For a generic potential the trajectories of test particles will not form closed non-circular 
orbits at all radii. Since we want to avoid intersections of the particle orbits, any deformation 
we make to the orbits in the disk must be a perturbation to a circular orbit. In this case we 
can use the epicyclic approximation, in which radial and vertical perturbations are harmonic 
oscillations with coordinate frequencies Uj. and ug, respectively. In this approximation, to 
first order in the perturbations e and 7, the map ([T]) can be written as 

(8a) 
(8b) 
(8c) 

The horizontal epicycles have amplitude e and phase x] the vertical epicycles have amplitude 
7 and phase (; in addition the azimuthal position is perturbed by 0o- The azimuth does not 
increase monotonically but undergoes a sinusoidal variation related to the radial oscillation 
through the conservation of angular momentum. In the Newtonian limit, r'^(f) = constant 
implies 

m = ^ . (9) 

For the Kerr metric in Boyer-Lindquist coordinates, conservation of angular momentum gives 

_ 2fu^{f) [1 + uj4f){a - sVGMf)] 
' ujr{f) f^-2GMf + a^ ' ^ ' 

where a is the black hole spin with a > for prograde orbits and a < for retrograde orbits. 

It's worth noting that in writing equations (IHl) we are neglecting time- independent per- 
turbations that arise from deforming one circular orbit into another. We could account for 
them simply by modifying the unperturbed disk. In our linearized treatment such pertur- 
bations are uninteresting because they cause no time variability. 
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In equations ([8]) we used the variable $(f, (p) = (f)(t = 0) = (j)+(j)o{f, (j))+P{f)e[f, $(f, x 
sinx[f, <l>(f, 0)], which can be solved iteratively for $, but such a solution is not necessary 
for our calculation. This expression can be rewritten using (I8bp as 

$(f, = (j) - uj^t - (3e[sm{uJrt + x) -sinx] • (U) 

The variables (f, $) give the original coordinates of particles in the perturbed disk. The 
choice of functional dependence of the perturbation functions on $ will allow us to write 
r(t) and 6{t) as functions of the mixed coordinates (f, 0) which will give the correct caustic 
structure provided that a certain initial condition is satisfied, as we will show below. 

The Eulerian coordinates r and 6 at fixed t must be periodic in 0. Therefore, using ([8]) 
we obtain 

e(f, $ + 2?™) = e(f, $), 7(f, $ + 2?™) = 7(f, $) , (12a) 

X(f, $ + 2™) = x(^, $) + 27rg(n), C(f, $ + 2?™) = C(f, $) + 27rg(n) , (12b) 

for any integer n. The function q{n) appearing in the two places above need not be one 
and the same, but we used the same symbol to economize notation. Since the perturbation 
functions must be smooth, this implies that q{n) depends only on n, and not on f or t for 
instance. 

In the epicyclic approximation ([8]), for any two functions M(f, $) and N{f, $) we obtain 

'dM\ fdN\ fdN\ fdM 



9$ 



dM\ fdN\ fdN\ fdM 



dr J ^\ J ^ \ 9r / ^ V 5$ 



(13) 



We do not want Jmn to become zero, hence we require {d^/d(l))f 7^ 0. This requires that the 
initial conditions in the perturbed disk are such that the angle $ = (f){t = 0) is a function of 
(j), otherwise from ([8]) we see that the map will no longer be one-to-one, which in turn would 
lead to caustics in the initial conditions. Since this term is time independent, once these 
initial conditions are satisfied, the zeros of Jmn and Jmn/ {d^/d(l))f will coincide. Thus, 
Jmn 7^ is a sufficient condition to justify our usage of (f, $) to label the particles in the 
disk. 

In complete analogy as above, we find that the ratio of Jmn to the Jacobian of the 
two functions with (f, 0) chosen as independent variables is equal to {d^/d(j))f,f From (ITT]) . 
{d^/d(f))f^t — 1 is first order (or smaller) in (3e and (5{de/ d^)f. Below we will show that for 
long-lived perturbations {de/d^)f = 0, and therefore the caustic structure is preserved if 
we restrict ourselves to small-amplitude perturbations (/3e ^ 1) and we choose to label the 
particles with coordinates (f, 0) or (f, $). 
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4. Perturbations in the Disk Plane 

Having developed the necessary background for kinematic density waves in possibly 
warped disks in flat or curved spacetime, we now ask whether there exist perturbations that 
last many orbital times before orbits intersect. We will find in general that orbits always 
intersect eventually in the absence of pressure or other non-gravitational forces. However, 
some patterns last many orbital times. We begin the examination in this section with 
perturbations that leave the disk fiat. 



4.1. Long-lived Patterns and their Lifetimes 

Perturbations in the plane of the disk satisfy T = Q = -n j1. Thus, from ([7]) we see 
that to avoid caustics we need J^p ^ 0. Let us write the Jacobian of the two functions 
[i?(f, 0; t), P(f, 0; if:)] in terms of the Jacobian of the functions [-R(f, $; t), P(f, $; t)] defined 
by composition with 0(f, $). With $ as the angular variable the calculation simplifies 
somewhat as the unknown function (j)o{f,(f)) is eliminated. 

Using (ISal) . (l8b|) and (IT5]) we can obtain Jjip/{d^/d(l>)f- The Jacobian has many terms, 
and to simplify things note that from f|T2l) . the derivatives of all perturbation functions with 
respect to f and $ are periodic in $. Also, we have to drop all terms proportional to e^, 
since our approximation does not conserve angular momentum at this order. However, this 
does not imply that we can drop terms such as e{de/df) and e{de/d^) since the derivatives 
need not be small. 

At t = the Jacobian evaluates to 

jRp{t = 0) ^ dx . de 

All partial derivatives with respect to f are taken holding $ constant. In order to avoid 
caustics in the perturbed disk initially we require Jrp ^ everywhere which will be satisfied 

*V+f^llV<l. (15) 



^dr ) \ dr 

Similarly, to avoid caustics forming within an orbit time it is sufficient to require 



If we assume that j3e e/r <^ 1 and the inequalities f|T5|) and f|T6|) are strongly satisfied. 
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then at late times the Jacobian becomes 



e ( - ^ J sm{uJrt + X) + -Q^ cos{uJrt + x) 



{d^/d(f))r 

where a prime denotes a derivative with respect to f, e.g. u'^ = du^/df, and 



(17) 



/i(f) ^ ^ . (18) 
To avoid intersection of orbits we require J/jp > 0. At late times this requires 

These conditions may be integrated to give x ~ /i(f)$ + 77(f) and e ~ e(f) where 77 is an 
integration constant. From (112bp . the first condition implies /x = m is an integer. Thus, 



to maximize the chance of having long-lasting waves we restrict perturbations to have the 
following form: 

e = e(f) , X = + v{^) ■ (20) 
The two functions e(f) and ri{f) specify the planar perturbations. 

From fl20l) we see that the first of conditions f|T9l) cannot be satisfied everywhere in the 
disk because /i is a function of f. Thus, intersections will inevitably occur: it is impossible 
to have an infinite-lived kinematic density wave except for special power-law potentials like 
Kepler and the simple harmonic oscillator. However, the timescale for caustic formation 
can be made arbitrarily large if one restricts perturbations to small neighborhoods around 
the discrete radii such that /i(rm) = m. Thus we consider a perturbation e(f) which is 
negligible everywhere except over a range 6r about r^. From (fT7|l and (l20l) . the perturbation 
lifetime until caustic formation is 

11 

In linear perturbation theory, 6r should be interpreted as 6f, i.e. the initial displacement of 
the unperturbed disk. However, it is more accurate to evaluate the orbital frequencies at the 
displaced position, in which case one can regard 6r as the maximum of the amplitude emax 
and the width 6f of the perturbation field e(r) (which might be, for example, a Gaussian 
centered at r^)- 

Thus, we have found relatively long-lived perturbations around discrete radii satis- 
fying fi{rm) = m = integer. Combining ([8]), (fTTI) . and (!20l) we find that to first order in e 
the perturbations have the form: 

r = f — e(f) cos[m0 + uJm{r)t + //(f)], Umir) = ^ri^) — ^^^{r) . (22) 
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This perturbation represents an |m|-armed spiral and Um is called the modulation speedj^ 
The modulation frequency of the spiral is fm = tUm/Svr, which is the frequency for an arm- 
to-arm displacement. Note that u'^ = at r = r.^; the modulation speed has an extremum 
at the radii of long-lived perturbations. 

In order to produce observable modulation, the pattern must not only survive against 
intersections, it must also not be too tightly wound. If, for fixed and t, r—f oscillates rapidly 
in the radial direction, the light curve will show little modulation. Significant modulation 
can persist only for a winding-up timescale defined by the condition that the argument of 
the cosine in fl22l) change by less than 2n as r changes for fixed (j) and t. At arbitrary f, this 
gives the condition + Sr) — co'm(f)] + ri'6f \ < 2ti. At f = r^, = the lifetime 

against winding is (if rj' is neglected) 

47r 27r 

The long lifetime has a simple interpretation. If the pattern speed ujp = —uj„i/m is constant 
over a range of radii, a density wave located over that range will remain constant in the frame 
rotating with angular velocity cUp, hence will avoid winding up. Note that if (5r ~ e (larger e 
is forbidden by the requirement e' < 1), the timescales for winding up and caustic formation 
are comparable; in general winding up precedes caustic formation. Had we retained nonzero 
r]', the winding up would have been faster on one side of r = r^, and slower on the other. 



4.2. Spiral Arm Direction 

In order to see whether the spiral arms are leading or trailing, we need to find the sign 
of (90/9r)spirai arms cvaluatcd aloug the spiral arms. Thus, we first find the initial spiral arm 
positions. These follow from the extrema of Jijp(f, $;t = 0). From (fT4|) and (!20|) . setting 
?7 = to maximize the perturbation lifetime, we obtain 

^ ^ I—- . (24) 

no 1 — e' cos(m$) 

Setting the gradient to zero gives conditions for extrema: e" cos(m$) = and me' sin(m$) = 
0. For generic e and m ^ this corresponds to the loci: (1) e' = and $ = (A^ -|- l/2)7r/m; 



2 



This is the angular speed of Ught curve oscillations, as opposed to the lower angular speed for the pattern 



to rot ate a full 2it. The corresponding pattern speed is tOp = —ujm/m = — cor/m. iBinnev fc Tremaine 
1987 ) define cOp — lo^ ~ {n/m)ujr where n and m are integers. This is exactly the result of lAmin fc Frolov 



20061 ). For n> 1 the orbits self-intersect so we do not consider this case. 
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and (2) e" = and $ = NTc/m; where is an integer. From the second partial derivatives 
we find that the first locus corresponds to saddle points in the density, while the second 
corresponds to overdensities if (— l)^e' > and — (— l)^e"' > at these points. If the signs 
are reversed, then we would have underdensities; and if the signs are mixed then we have a 
saddle. An example of this result is given in the first snapshot of Figure [3] for the m = —2 
mode in the Kerr metric. There we have a gaussian e(f) for which we obtain four minima 
and four maxima in the density. If at each radius we maximize n/riQ with respect to $, we 
find that initially the spiral arms are radial segments lying at $ = A^vr/m. Thus, there are 
2\m\ overdensities for the mode with azimuthal number m. 

The pattern makes a 2tt revolution with a period 2TTm/ujjn, which implies that (p + 
Umt/^ = Nn/m along the spiral arms. Since we need only the sign of {d(f)/dr)spira.i arms, we 
can work to zeroth order in e. The sense of winding up of each spiral arm is then given by 

— =uJpt = = ^(^-m). (25) 

Thus, assuming that co'm(^) is monotonic across the perturbation, we can conclude that 
overdensities which he inside r.m have opposite sense of winding than those lying outside. 
This is consistent with our earlier result that the pattern speed changes sign across the radii 
Tm- This behavior is exemplified Figure [3] for the m = —2 mode in the Kerr metric. Leading 
waves correspond to a pattern speed increasing with radius, or u'^ < 0. 



4.3. Spherical potentials in the Newtonian limit 

Applying the results obtained above requires finding integer solutions of /i(f) = m. 
In the Newtonian limit, with gravitational potential V{f), the effective potential for radial 
motion is V^fr(^) = V{f) + L^/2f^. For a unit mass test particle we can then write 



L 



^V'(f), = 7(K'4(^1)L=const = ^j^y'ir) + V"{f) 



1^ = — 



-V'(f 

r 



V'{f) -fV"{f)_ 



(26) 
(27) 



in the Newtonian limit. Note that the derivatives in ix = ui'^/ui'^ are taken across particle 
orbits and thus are not evaluated at L =const. 



As an example, we consider a set of classical potentials used in the literature. Surpris- 
ingly, we find that they are divided into three classes: 
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1. Potentials which do not have integer solutions for /i. These include: power law poten- 
tials V{r) oc — 1/r" with 2 > n > 1, the Hernquist potential V{r) = —a/{r + b), the 
Jaffe potential V{r) oc ln(r + a), and the NFW potential V{r) = —{a/r) ln(l + r/b), 
where a and b are constants. These potentials admit no long-lived kinematic density 
waves. 

2. Potentials with one solution for /i = 2 (and another one for /i = 3 at r = 0): the 
potential V{r) = —{a/r)\n[r/b + ^/ q + {f/by], the logarithmic halo potential V{r) oc 
ln(r^ _l_ q2-j^ ^Yie Plummer potential V{r) oc — l/Vr^ + a^. 

3. Potentials with non-positive integer solutions for fi: V{r) = —b[l + (a/r)^]/r, V{r) oc 
1 — exp(a/r), and the effective potential of the Kerr metric, which is discussed below. 



5. Vertical and Three-dimensional Perturbations 

Next we consider perturbations that warp the disk out of its plane, starting with the 
case e = but 7 7^ 0. The discussion of vertical perturbations follows closely the earlier 
presentation. The relevant Jacobians are given by: 

■Jrt dC, 

— cos {uet + C) - iT^ sin {uet + C) , (28a) 



J BP 



{d^/d(t))r 

JpT 

{d^/d^)f 



1 (28b) 
£ cos {ujet + C) + sin {ujet + () , (28c) 



where 



For orbits in a spherical potential, uJe = uj^ implying u = 1 for all r. For a nonspherical 
potential, z/ 7^ 1. 

As long as we require that (5$/90)f 7^ 0, we have Jrp 7^ 0, implying that no caustics 
will ever form since at least one contribution to ([7]) is nonzero. By comparing (fT7|) with 
f l28cp . we can expect to find a spiral pattern for the vertical perturbations. Let us find and 
maximize the winding-up timescale. This occurs when the oscillations of Jpt are bounded 
in amplitude as t — 00. By analogy with equations (HM and (12(1 we obtain 



7 = lir) , C = n$ + ^(r) , 



(30) 
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where n is an integer. The two functions 7(r) and S,{r) specify the vertical perturbations. 
Similarly to equation fl22l) . we find 

9 = 9- 7(f) cos[n0 + Qrn{r)t + ^(f)], n„(f) = uje{f) - nuj^{f) . (31) 

As before, the pattern is an |?T,|-armed spiral whose modulation speed has an extremum at 
the radii of long-lived perturbations: VL'^ = at r = r„. 

The vertical winding-up timescale is defined by the condition that the argument of the 
cosine in (1311) changes by less than 2tt as r changes for fixed and t. By analogy with the 
planar case, near the critical radii where v{f) = n we find (for nonspherical potentials and 

e' = o) 

An 271 

^ KiW ^ ^yW? ■ ^^^^ 

A nonzero ^' would decrease on one side of r„. For a spherical potential, Vtm = for 
n = 1 corresponding to a stationary warp. Such a perturbation is long-lived but provides 
no modulation for a light curve. A weak nonspherical distortion of the potential or a weak 
non-gravitational torque would produce a nonzero pattern speed, hence modulation, and 
might still preserve a long lifetime against winding up. 

Note that the Jacobian Jjix (I28al) is well-behaved and does not lead to further con- 
straints. 

Now let us consider general orbit perturbations with 7 7^ 7^ e. The first question that 
arises is whether we can still use the previous results. To see that this is indeed the case, let 
us find the relevant three Jacobians. The Jacobian Jp^ is exactly the same as in (I28cl) . and 
the Jacobian J^p is the same as the one in f[T7|) . The Jacobian Jrt is quite complicated in 
general, but it simplifies in the limit of large times: 

'^^{di/d(^^ = eju'^tifin - um) sin(cj^t + x) sm{uj0t + () , (33) 

where /i = oj'^juj'^ and v = oj'q/oj'^ as before. The integers m and n follow as before from 
periodicity conditions on x ^^id C,. 

The condition for caustics to form is that all three Jacobians vanish at the same point. 
In general this is impossible. However, perturbations can become tightly wound, after which 
time they would produce little modulation in a light curve. Maximizing the winding up time 
leads to the epicyclic solutions ( !22l) and ( !3T|) found before. 

A necessary (but not sufficient) condition for the three-dimensional wave to be long- 
lived (i.e., the inverse timescale is quadratic in 5r /r) is that the perturbations are centered 
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around radii Vc which satisfy both /i(rc) = m and z/(rc) = n. The requirement that both 
/i(rc) and //(r^) are integers is very restrictive and is possible only for spherical potentials 
where u = 1 for all r, or other special cases. 



Accretion Disks Around Kerr Black Holes 



The Kerr metric in Boyer-Lindquist coordinates is stationary and axisymmetric with 
orthogonal metric on constant time slices: 

(34) 



ds"^ = -a^dt^ + w'^{d(t) - ujdtf + ^rfr^ + p'^dO'^ . 



A 

In geometrized units G = c = 1, for a black hole of mass M and spin a (with —M < a < M 
to allow for prograde or retrograde orbits), 

2Mr 



=r'^ + cos^ 6* , A = - 2Mr + , ~l 



2Mra 



p2A + 2Mr(a2 + r2) 



zu 



? 



sin^ e . 



(35) 



Throughout this section, when we give frequencies or timescales in SI units, those are given 
for a one solar mass black hole. The physical frequency is obtained by multiplying the value 
for the solar-mass black hole by Mq/M. 

On shces of constant coordinate time {dt = 0), the functions A, B, and C in ([3]) can be 
read off from the metric: 



A 



P 



, B = w, C = p . 



(36) 



As desired, these quantities are nonzero outside the event horizon except at the poles (a coor- 
dinate singularity). The epicyclic coor dinate frequencies for circular orbits in the equatorial 
plane of the Kerr metric are given by (jPerez et al.l 119971 ) 

1 



M [M 
[1- 6 h Sa\ — - 3— 



1/2 



7^ V .y^2 



[M 3a2 
u;e = iJ^ I l-4aW— + — 



1/2 



(37a) 
(37b) 
(37c) 
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Using fl371) we can compute both fi and z/. The equation fi{rm) = m has integer solutions 
for all m < 0, with m ^ 1^ as r —>■ oo, recovering the result for the Newtonian Kepler 
potential. From fl22]) we conclude that the spiral patterns formed at move in the prograde 
direction (i.e., in the same direction as the particle bulk flow) for m < 0. The radii are 
plotted versus the black hole spin parameter, a, in Figure [H Notice that approaches the 
innermost stable circular orbit (ISCO) for large — m. The properties of the m = mode 
are qualitatively different from those of the modes with m ^ 0, and this can be seen in all 
plots that follow. The reason is that when m = the pattern frequency is determined by 
the radial oscillation frequency, while for the rest of the modes Um is determined mainly by 
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Fig. 1. — The radii for planar density waves satisfying fi{rm) = "rn, versus black hole 
spin for the Kerr solution. The integer mode number m is zero for the uppermost curve and 
then decreases with decreasing r^. As m ^ — oo, approaches riscO) the innermost stable 
circular orbit. 




Fig. 2. — The modulation frequency (the frequency of observed light curve variations) for 
planar density waves of different m in the Kerr metric. The lowest frequency curve corre- 
sponds to m = 0, with the frequency increasing with |m|. These frequencies appear to be 
too large to account for observed black hole QPOs. 
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Figure [2] shows the modulation frequency fm = i^m/27r. The lowest mode m = has 
a frequency of fm > 700 Hz, and the next mode has a frequency exceeding 2 kHz. One 
can com pare these frequencies with mea sured QPO frequencies for well-studied black hole 



binaries (IRemillard &: McClintockl |2006| ). The lower of the two commensurate frequencies 



are 300 Hz for GRO J1655-40, 184 Hz for XTE J1550-564, and 113 Hz for GRS 1915+105; 



for ea ch case there is a second frequency about 1.5 times higher. IRemillard fc McClintock 



( I2OO6I ) show that the lower of the two frequencies is approximately 1860 Mq/M Hz for these 
systems, too small to be consistent with the modulation frequencies of Figure [2] for |m| > 0. 
This is a shame, given that the frequency ratios for m = —2 and m = —3 have nearly a 3:2 
ratio over a wide range of black hole spin. Thus, the simplest model of kinematic density 
waves fails to account for the black hole high-frequency QPOs. 

Our results for the pattern speed and of the m-th mode exactly reprod uce the fre- 



uency and radius of the same g-mode as derived in the diskoseismology model of iPerez et al. 



19971 ) in the pressureless limit (cf. our Figure [2] and their Figure 5). In the limit of large 
vertical mode number, the frequency difference between our and the diskoseismic result is 
proportional to the sound speed and the radi al mode nun i ber, a nd is inversely proportional 



to the vertical mode number [as indicated in iPerez et al.l (Il997l ). equation (5.3)]. Since we 



consider only pressureless disks, we recover no analogs of p- and c-modes. However, our 
analysis of g-modes in pressureless disks is valid for generic potentials, and it allows for 
the interpretation of the perturbations as kinematic density spiral arms, as well as for the 
visualization of the patterns that form in the disk. Thus, we are able to see that the m-th 
mode has 2\m\ overdensities (see below) resulting from pure kinematics, which might dou- 
ble the expected frequency of the X-ray modulations. This result has been mi ssed by the 



diskoseismic models. The reason for that is that in the diskoseismic analysis, iPerez et al. 



( I1997I ) assume a time dependence of the Lagrangian perturbations of the form e**^*, where a 
is independent of f. Comparing this with our equation fl8al) . we can read off the Xds for the 
diskoseismic models: 

Xds = m$ + t{a - ujr{f)) + /3eO(l). 

Thus, the choice of time dependence in the diskoseismic analysis breaks the epicyclic ap- 
proximation after a timescale of 



which is ~ (5r)~^ times shorter than the lifetime of the perturbations that we find. In the 
pressureless limit, nothing can break the validity of the epicyclic approximation. In the 
diskoseismic models this is reflected in the fact that the radial width of the perturbations 
vanishes in that limit. Therefore, our result generalizes the result of the diskoseismic models 
for pressureless disks by allowing for transients and maximizing their lifetimes. 



Fig. 3. — Spiral-arm pattern caused by a perturbation with m = —2 around a Kerr black hole 
with a spin parameter a = 0.5M. Time advances from left to right and top to bottom. The 
dark circle shows the circle where fi{r) = —2, i.e. r_2 = 4.33M. The ISCO is risco = 4.23M. 
The spiral pattern is formed by a small-amplitude perturbation in the plane, as described 
in the text. Since the perturbation is nonzero only in an extremely thin annulus around the 
black hole, we have cut out the central portions of the disk for clarity. Thus, the center of the 
polar plot does not correspond to r = 0, but to r = r_2 — 0.015M, and the plot extends to 
r = r_2 -I- O.OIM. The last two snapshots already have formed caustics. Note that although 
m = —2, there are in fact 4 overdense regions because the waves are leading inside r_2 and 
trailing outside r_2. 
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To illustrate the kinematics of the planar density waves, in Figure [3] we plot the m = —2 
spiral pattern for several times. The initial perturbation was chosen with //(r) = and 
e(r) = 10~^Mexp[— (r — r_2)^/2cr^], with (5r ~ a = 0.003M. The arm-to-arm displacement 
period is 0.87 ms. The timestep between the snapshots is 9.3 s, and we have made it to 
be an exact multiple of the pattern frequency. Therefore, the pattern appears stationary, 
although it rotates counterclockwise. The winding-up timescale is about 30 s, and the 
caustic-formation timescale is about 1 minute. We ran similar calculations for a set of 
perturbations with different m, a, e and t], and we found that our analytical estimates for 
the winding-up timescale and the lifetime of the pattern are consistent with the numerical 
results. 

Figure [3] verifies the analytic expectations of section 4.2. We see that for rj = the spiral 
pattern for a mode with mode number m has 2\m\ overdense regions — \m\ trailing waves 
on the outer edge of the perturbation and |m| leading waves on the inner edge in agreement 
with equation (l25l) . If these overdensities produce comparable modulation of the intensity, 
then an observer would see a modulation at a frequency 2/^, exacerbating the problem of 
matching kinematic density waves to observed black hole QPOs. 

If simple density waves cannot explain the observed QPOs, maybe a beat frequency 
between an orbiting blob and density waves can do the job. Consider a blob of gas orbiting 
the black hole near with angular frequency uj^{rjn)- It must be small and dense enough 
not to be tidally sheared apart. For m < —1 the spiral arms will rotate faster than the blob 
around the disk and thus the blob will be compressed every time the spiral waves pass over 
it. This leads to a beat frequency in the emission from the blob given by {um — \m\uj,j))/2'n' = 
u!r/2TT = fr{rm), i.e. the radial epicyclic frequency. While the blob remains on one side of 
Tm, its beats lead to a flickering frequency of fr- If the blob straddles hence is overtaken 
by both leading and trailing density waves, the flickering frequency is 2/^. 
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Fig. 4. — The radial oscillation frequency at versus the black hole spin parameter. The 
uppermost curve corresponds to m = 0, with decreasing as m becomes more negative. 
The beat frequency between a blob in the disk and a kinematic density wave near is 
either or 2/^ (the frequency being doubled if the blob straddles r^). These frequencies 
appear to be too small to account for observed black hole QPOs. 
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The radial frequency fr is plotted in Figure |H The radial frequency in the Kerr metric 
is less than the azimuthal frequency and goes to zero at the ISCO (corresponding to modes 
with m —>■ — oo). This frequency is too low to account for high-frequency QPOs, even if it is 
doubled. 
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a/M 

Fig. 5. — Low-m sets of frequencies having approximately a 3 : 2 ratio. The curves are 
labeled by the appropriate ratios of /m,|m| for density waves or fr,\m\ foi' blob-wave beat 
frequencies (in all cases \m\ < 0). Ratios of fm/ fr exceed 2.0 except for the special case 
/m,o = /r,o- Although the pairs of modes with m = — 2 and m = — 3 have ratios close to 1.5 
over a wide range of spin parameter, the actual frequencies are discrepant with observations. 
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Maybe the black hole mass estimates are in error. If so, the predicted frequencies are 
wrong but frequency ratios are robust. Thus in Figure [5] we plot all frequency ratios of 
low-m modes that are close to 1.5. Both the density wave modulation frequencies fm and 
the radial (blob-wave beat) frequencies fr are included. The only ratios that come close to 
1.5 are the pairs with m = —3 and m = —2. However, the density wave frequencies are too 
high (the lower frequency is at least 2500 Hz for one solar mass) and the beat frequencies 
are too low (the lower frequency is less than 420 Hz or 840 Hz if it is doubled) compared 
with the observed range 1100 to 2000 Hz. The nearest frequency match occurs for spiral 
density waves with m = —2 for a/M = 0, which would require that the black hole masses be 
about 25% larger than the upper range of observational estimates. For GRS 1915+105, for 
example, the lower 113 Hz QPO could be explained by fm,~2 if the black hole mass were 22 
Mq. The frequency ratio 168/113 is close to /m,-3//m,-2 = 1-47. However, the spiral density 
wave model also predicts modulation at frequencies twice as high (from the combination of 
inner leading and outer trailing arms), which has not been seen. 
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Fig. 6. — The winding up timescale for spiral density waves (the uppermost curve corresponds 
to m = 0), giving an order of magnitude estimate of their hfetime in minutes. The hfetimes 
are much longer than they need to be to produce the observed widths of the QPO spectral 
peaks. 
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For completeness we present in Figure [6] the lifetimes of various modes to winding up. 
The lifetimes must be at least 10 times the inverse frequency in order that mode decay not 
broaden the peaks in the power spectrum to a lower Q value than is observed. However, 
this condition is very easily satisfied for modes concentrated near the discrete radii r^- The 
broadening of the timing peaks is probably not due to finite mode lifetime unless the modes 
extend in radius a significant fraction of the distance to the next discrete mode radius r^. 

In addition to the beating of a blob with a wave, one might consider the nonlinear 
interaction between waves of different pattern speeds. However, this leads nowhere because 
the pattern (or modulation) speed is a function of radius, and spiral density waves of different 
mode number have different frequencies only because they are localized at different radii. 

The application to Kerr has until now considered only planar perturbations. Vertical 
perturbations can also be excited and are long-lived when concentrated at the discrete radii 
Tn such that z/(r„) = n. For the Kerr metric there are only two solutions outside the ISCO: 

1. For a/M > 0.9524 there exists a solution with = at tq ~ 1.9M. The frequencies 
are fg = Qm/^n = 4.9 kHz, with a winding-up timescale of about r = 2 min for 
6r ~ O.OIM. 

2. For a/M > 0.9987 there exists a second solution with n = —1 at r_i ^ 1.2M. The 
frequencies are fg = 2.7 kHz, Qrn/2n = 16.5 kHz, with a winding-up timescale of about 
T = 0.4 min for 5r ~ O.OIM. 

Because these modes occur outside the ISCO only for a/M close to maximal, the mode 
radii and frequencies are nearly constant over the small range of black hole spin where they 
exist. Given the very narrow range of spins for which these long-lived warps exist, and given 
their very high frequencies, they are likely to be observationally unimportant. 

At larger radii and for smaller black hole spins, the lifetime of vertical waves becomes 
relatively large even though fi and u are not integers. To order of magnitude the winding time 
multiplied by mode frequency is ~ fiuJmTw ~ {r/6r) and similarly for the vertical modes. 
Also, the horizontal and modulation frequencies become comparable: for the same mode 
number, 0.9 < Um/^m < 1 for radii corresponding to pattern frequencies in the range 10~^ 
Hz to 10^ Hz (multiplied by Mq/M). This frequency ra nge includes the low- frequency QPOs 



seen in the soft power-law state of black hole binaries (IRemillard fc McClintockll2006l ). The 
QPO frequency varies with X-ray flux, implying that the modulation does not arise at a 
fixed radius (unlike the high-frequency QPOs). 

One explanation for the low-fequency QPOs is Lense-Thirring precession (Stella et al. 
1999, Schnittman et al. 2006), with modulation speed VL^t = ^<i> — ^e- At large radii, z/(r) 
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1, so vertical density waves with n = 1 are long-lived and have pattern speed equal to ^lt- 
However, planar perturbations have /i(r) 1 and are comparably long-lived. Therefore, 
either planar (spiral waves) or vertical waves (warps) are candidate explanations for low- 
frequency QPOs, and they have practically the same frequencies and lifetimes. Instead 
of having no density wave solution consistent with observations (the situation for high- 
frequency QPOs), the low- frequency QPOs have a plethora of kinematic solutions with little 
to distinguish them. 



7. Conclusions 

Starting with a simple kinematical consideration — fluid streamlines should not intersect 
— we were able to classify all long-lived linear perturbations in pressureless thin accretion 
disks in generic potentials including the Kerr metric. The result is a set of patterns (planar 
spiral arms, vertical warps, and linear combinations) whose mode functions are concentrated 
at discrete radii satisfying a simple condition: the pattern frequency for non-self-intersecting 
orbits must be approximately constant over a range of radii, i.e. its derivative with respect 
to radius must vanish. This condition automatically picks out discrete set of frequencies 
for modes with a long lifetime against winding up. The lifetimes are 0(r/5r)^ oscillation 
periods where 6r is the radial extent of the mode. 

Although our treatment neglected p ressure, the r esults are in excellent agreement with 



the diskoseismic (7-modes investigated by lPerez et al.l (Il997l ). when the sound speed is much 



less than the orbital speed. Our treatment, based on the epicyclic approximation, is much 
simpler than the solution of the relativistic fluid equations. 

One of our goals was to investigate the possibility of explaining the 3 : 2 ratio in high- 
frequency quasi-periodic oscillations observed in accreting black hole binaries. Our model 
gives several mode pairs with a frequency ratio close to 1.5 for a wide range of black hole spin. 
This result is a direct consequence of the discrete radii at which the spiral patterns exist. 
However, we cannot recover naturally the correct frequencies unless the black hole mass 
estimates are seriously in error. Planar spiral density waves have frequencies too high, while 
the beat frequency between a blob and these waves is too low, compared with the observed 
frequencies. Moreover, the best matches follow for modes with |m| = 2 and \m\ = 3, with 
no explanation for the failure to observe the mode with |m| = 1. 

Because our analysis classified all long-lived perturbations of pressureless disks, the fail- 
ure to reproduce the observed high-frequency quasi-periodic oscillations implies that pressure 
or other non-gravitational forces (e.g. magnetic stresses or radiation pressure) must be im- 
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portant in producing the QPOs. 

Based on the diskoseismic results, we would expect that small non-gravitational forces do 
not significantly change the mode frequencies. Large non-gravitational forces would thicken 
the disk and might destroy the coherence and longevity of modes, although this remains 
to be seen. Even so, some mechanism is needed to pick out discrete frequencies of the 
high-frequency QPOs. 

These results make more attractive the possibility that nonlinear coupling between 
modes, e.g. parametric resonance, is responsible for selecting the frequency ratios (Abramow- 
icz & Kluzniak 2001; Rebusco 2004). One approach to investigating this process is to carry 
the perturbation methods of the current paper to one higher order, to include a coupling 
between radial and vertical oscillations. 

This work was supported by NASA grant NAGS- 13306. 
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